/*****************************************************
 * 这是一个矩阵计算库，由于自己之前用eigen库非常不方便，我参考了翀哥代码也复现了一下矩阵计算库
 * 作者：陈冠宇
 * 特别鸣谢：黄文翀
 *****************************************************/

#ifndef MATRIX_H
#define MATRIX_H

#include <iostream>
#include <cmath>
#include <cstring>

const double mat_eps = 1e-12;

class Matrix;
class ColVector;
class RowVector;

class Matrix
{

private:
    double *a;

public:
    int n, m;
    Matrix();
    Matrix(const int &_n);
    Matrix(const int &_n, const int &_m);
    Matrix(const Matrix &A);
    Matrix(const double *p, const int &_n);
    Matrix(const double *p, const int &_n, const int &_m);
    ~Matrix();
    inline bool is_empty() const { return n == 0 || m == 0; }
    inline bool is_square() const { return n == m; }
    Matrix &operator=(const Matrix &rhs);

    const double element(const int &r, const int &c) const; // 用于方便调用数组元素
    double &element(const int &r, const int &c);
    double vecnorm(const double &p) const;
    const double operator()(const int &r, const int &c) const;
    double &operator()(const int &r, const int &c);

    RowVector getRow(const int &r) const;
    ColVector getCol(const int &c) const;

    Matrix operator+(const Matrix &B) const;
    Matrix operator-() const;
    Matrix operator-(const Matrix &B) const;
    Matrix operator*(const Matrix &B) const;
    Matrix operator/(const double &p) const;
    Matrix T() const;
    void swaprow(const int &r1, const int &r2);
    void swapcol(const int &c1, const int &c2);

    friend Matrix diag(const Matrix &A);
    Matrix getSubmatrix(const int &r1, const int &r2, const int &c1, const int &c2) const;
    void setSubmatrix(const int &r, const int &c, const Matrix &rhs);

    friend ColVector solve(Matrix A, ColVector b);
    double det() const;
    Matrix inv() const;
};
/***************************列向量****************************/
class ColVector : public Matrix
{
public:
    ColVector() : Matrix() {}
    ColVector(const int &n) : Matrix(n, 1) {}
    ColVector(const double *p, const int &n) : Matrix(p, n, 1) {}
    ColVector(const Matrix &rhs);
    inline int size() const { return n; }
    const double operator()(const int &inx) const;
    double &operator()(const int &inx);
    ColVector &operator+(const ColVector &rhs) const;
    ColVector &operator-(const ColVector &rhs) const;
    ColVector &operator-() const;
    RowVector T() const;
};
/***************************行向量****************************/
class RowVector : public Matrix
{
public:
    RowVector() : Matrix() {}
    RowVector(const int &m) : Matrix(1, m) {}
    RowVector(const double *p, const int &n) : Matrix(p, 1, n) {}
    RowVector(const Matrix &rhs);
    inline int size() const { return m; }
    const double operator()(const int &inx) const;
    double &operator()(const int &x);
    RowVector &operator+(const RowVector &rhs) const;
    RowVector &operator-(const RowVector &rhs) const;
    RowVector &operator-() const;
    ColVector T() const;
};

/***************************线性方程组计算（数值代数）****************************/
double det(const Matrix &A);
Matrix inv(const Matrix &A);
ColVector solveLowerTriangular(const Matrix &A, const ColVector &b);
ColVector solveUpperTriangular(const Matrix &A, const ColVector &b);
Matrix choleskyImproved(const Matrix &A);
ColVector solveByLDL(const Matrix &A, const Matrix &b);

/******************************相关运算符重载******************************/
Matrix operator*(const double &k, const Matrix &A);
ColVector operator*(const double &k, const ColVector &x);
RowVector operator*(const double &k, const RowVector &x);
ColVector operator*(const Matrix &A, const ColVector &x);
RowVector operator*(const RowVector &x, const Matrix &A);
double operator*(const RowVector &rv, const ColVector &cv);
std::istream &operator>>(std::istream &in, Matrix &A);
std::ostream &operator<<(std::ostream &out, const Matrix &A);

Matrix zeros(const int &n, const int &m);
Matrix eye(const int &n);
ColVector zeroCol(const int &n);
RowVector zeroRow(const int &n);

#endif
